/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined _weight_set_operator_
#define _weight_set_operator_

#include <src/data/ParticleSet.h>
#include <src/miscellaneous/OS_basic.h>
#include <src/data/Matrix.h>
#include <src/math/Operators/GenericAlgebra.h>
#include <src/test/miscellaneous/stat/data_performance_analyzer.h>

/*! \brief default interaction for class weighted_operator
*/
template <class A, std::size_t Dim = A::Dim>
struct No_Interaction {
	typedef typename A::value_type value_type;
	typedef typename A::const_reference const_reference;
	typedef typename A::reference reference;
	//! fill an already initialized reference with \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. 
	void operator()(const_reference left, const_reference right, reference out) const
	{
		out = right;
	};
	//! return \f$ \frac{S_i-S_j}{|S_i-S_j|^2} \f$. slow method, use 
	value_type operator ()(const_reference left, const_reference right) const
	{
		return right;
	};
};

template <class T, class I>
struct weighted_operator_interface 
	: public virtual graph_base<typename T::argument_base_type>
{
	typedef I interaction_type;
	weighted_operator_interface(const weight_set_type & from) : graph_base(from) {
	};
	weighted_operator_interface() {};
	virtual ~weighted_operator_interface(){};

	virtual typename T::result_type operator()(typename T::argument_type from) const 
	{
		typedef typename T::argument_base_type argument_base_type;
		typedef typename T::argument_type argument_type; 
		typedef typename T::result_type::value_type value_type; 
		typedef typename argument_base_type::reference reference; 
		typedef typename T::argument_base_type argument_base_type; 
		typedef typename T::result_type result_type; 
		typedef argument_base_type::size_type size_type;
		typedef argument_base_type::const_reference reference_argument_type;
		typedef result_type::reference reference_result_type;
		typedef argument_base_type::const_iterator const_iterator;
		typedef result_type::iterator iterator;

		size_type i,j;
		result_type result(from.size_list());
		result.assign(0);

		const_iterator beg = from.begin();
		value_type temporary(*from.begin());

		typedef weight_set_type::const_iterator1 const_iterator1;
		typedef weight_set_type::const_iterator2 const_iterator2;
		const_iterator1 itl;
		const_iterator2 itl_;
		const_iterator1 end(weight_set_.end1());
		for (itl = weight_set_.begin1();itl != end;++itl) {
				i = itl.index1();
				reference_argument_type ref_i(*(beg+i));
				const_iterator2 end2(itl.end());
				for (itl_ = itl.begin();itl_ != end2;++itl_) {
						j = itl_.index2();
						reference_argument_type ref_j(*(beg+j));
						interaction_(ref_i,ref_j, (reference) temporary);
						temporary *= (*itl_);
						result[i] += temporary;
				};
				
		};
		return result;
	};
	virtual void output (std::ostream& o) const {
		o<<weight_set_ << std::endl;
	};
	interaction_type interaction_;
};


/* /brief Weighted Operators algebra

	The class weighted_operator implement an algebra of Weighted Operators (inner (+,-,*), outer (*,/)).
	We call a Weighted Operators a function f that act as \f$ R = f(S) \f$, where
	\f$ R=\{R_i\}_i \f$ and \f$ S=\{S_i\}_i \f$ are collections, and f acts as
	\f[
		R_i = \sum_{j \in \nu_i} \omega_{i,j}I(S_i,S_j)
	\f].
	- The set A is described by the template parameter A.
	- The set R is described by the template parameter R.
	- The set \f$ \omega_{i,j}_{\{i,j\}} \f$ is of kind D. We use 
	a sparse matrix container for this set. The default implementation is for the default container weight_set_default_def.
	- The interaction function I is described by the template parameter I. Default behavior is no interaction at all (finite difference operator):
	\f[
		R_i = \sum_{j \in \nu_i} \omega_{i,j} S_j
	\f].
	
*/
template <class A,class R, class I = No_Interaction<A>, class V = smart_ptr<A> >
struct weighted_operator
	: public OS_Generic_Algebra<A , R, V> 
	,	boost::additive< weighted_operator<A,R,I,V> 
	,	boost::multiplicative< weighted_operator<A,R,I,V>, typename outer_structure<A,R>::base_value_type
	> >
	, public weighted_operator_interface< OS_Generic_Algebra<A , R, V>, I >
{
	enum{Dim = A::Dim};
	typedef weighted_operator									type;
	typedef typename A::Boundary_type							Boundary_type;
	typedef OS_Generic_Algebra<A, R, V>							Base;
	typedef typename Base::argument_base_type					argument_base_type;
	typedef Base												Base_Cast;
	typedef weighted_operator_interface< Base, I >				interface_type;

	weighted_operator(){};
	weighted_operator(const weight_set_type & data) : graph_base(data){};
	weighted_operator(const Base &first) : Base(first)
	{};
	weighted_operator(const weighted_operator & from) : graph_base(from.weight_set_), Base(from) {
	};
	weighted_operator & operator = (const weighted_operator & from) {
		weight_set_ = from.weight_set_;
		Base::operator =(from);
		return *this;
	};
	virtual ~weighted_operator(){};

	type & operator-(){
		Base_Cast::operator-();
		weight_set_ *= -1.;
		return *this;
	};
	type & operator+=(const weighted_operator & from ){
		Base_Cast::operator+=(from);
		weight_set_ += from.weight_set_;
		return *this;
	};
	type & operator-=(const weighted_operator & from ){
		Base_Cast::operator-=(from);
		weight_set_ -= from.weight_set_;
		return *this;
	};
	type & operator*=(const weighted_operator & from ){
		Base_Cast::operator*=(from);
		weight_set_ *= from.weight_set_;
		return *this;
	};
	type & operator *=(const base_value_type & from) {
		Base_Cast::operator*=(from);
		weight_set_ *= from;
		return *this;
	};
	type & operator /=(const base_value_type & from) {
		Base_Cast::operator/=(from);
		weight_set_ /= from;
		return *this;
	};
	virtual result_type operator()(typename argument_type from) const 
	{
		if (first_)
			return (*first_)(from);
		else
			return interface_type::operator()(from);
	};

	virtual smart_ptr_type Create_Instance() const {
		return Create_Instance_helper_arg<type,weight_set_type>(&weight_set_);
	};
	virtual bool equals (const Base& from) const 
	{
		if (weight_set_ != static_cast<const type&> (from).weight_set_)
			return false;
		return Base::equals(from);
	};
	virtual void output (std::ostream& o) const {
		interface_type::output(o);
		Base::output(o);
	};
};

/* /brief Identity for weighted Operators 

	This class is the identity for weighted operator  :
	\f$ R=\{R_i\}_i \f$ and \f$ S=\{S_i\}_i \f$ are collections, and f acts as
	\f[
		 \{f(S_i)\}_i = \{S_i\}_i
	\f].
*/
template <class A>
struct weighted_operator_identity
	: public weighted_operator<A,A> 
{
	typedef weighted_operator_identity								type;
	typedef weighted_operator<A , A>								Base;
	typedef typename Base::Base_Cast								Base_Cast;
	typedef weighted_operator_identity								type;
	typedef std::size_t												size_type;		
	weighted_operator_identity(size_type size) : size_(size){
		weight_set_.resize(size_,size_,false);
		for (size_type i = 0; i < size_;++i)
			weight_set_(i,i) = 1;
	};
	weighted_operator_identity(const weighted_operator_identity &from) : size_(from.size), Base(from)
	{
	};
	weighted_operator_identity(const Base & from) :  Base(from) {
		OS_DYNAMIC_CHECK(dynamic_cast<const type&> (from) != 0," weighted_operator_identity::weighted_operator_identity : is this normal ?");
		size_ = static_cast<const type&> (from).size_; 
	};
	virtual ~weighted_operator_identity(){};
	virtual typename Base_Cast::smart_ptr_type Create_Instance() const {
		return Create_Instance_helper_arg<type,size_type>(&size_);
	};
	virtual bool equals (const Base_Cast& from) const 
	{
		if (size_ != static_cast<const type&> (from).size_)
			return false;
		return Base_Cast::equals(from);
	};

	Base & operator-(){
		Base_Cast::operator-(*this);
		return *this;
	};
	Base & operator+=(const weighted_operator & from ){
		Base_Cast::operator+(*this,from);
		return *this;
	};
	Base & operator-=(const weighted_operator & from ){
		Base_Cast::operator-(*this,from);
		return *this;
	};
	Base & operator*=(const weighted_operator & from ){
		Base_Cast::operator*(*this,from);
		weight_set_ *= from.weight_set_;
		return *this;
	};
	virtual typename Base_Cast::result_type operator()(typename Base_Cast::argument_type from) const 
	{
		return from;
	};
	std::size_t size_;
};



#endif